Approximation  Methods  for  Inverse  Problems  Governed 
by  Nonlinear  Parabolic  Systems* 

H.  T.  Banks]  C.  J.  Musante]  and  J.  K.  Raye 
Center  for  Research  in  Scientific  Computation, 

and 

Department  of  Mathematics 
North  Carolina  State  University 
Raleigh,  NC  27695-8205 

December  17,  1999 


Abstract 

We  present  a  rigorous  theoretical  framework  for  approximation  of  nonlinear  parabolic 
systems  with  delays  in  the  context  of  inverse  least  squares  problems.  Convergence 
of  approximate  optimal  parameters  and  that  of  forward  solutions  in  the  context  of 
semidiscrete  Galerkin  schemes  are  given.  Sample  numerical  results  demonstrating  the 
convergence  are  given  for  a  model  of  dioxin  uptake  and  elimination  in  a  distributed 
liver  model  that  is  a  special  case  of  the  general  theoretical  framework. 
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1  Introduction 


In  this  paper  we  present  results  for  approximation  of  parameter  estimation  problems  governed 
by  nonlinear  parabolic  partial  differential  equations  with  delays.  Motivated  by  a  concrete 
example  for  dioxin  uptake  and  transport  in  a  spatially  distributed  model  of  the  liver,  we  for¬ 
mulate  an  inverse  problem  in  an  operator  theoretic  setting  for  a  least  squares  optimization 
problem.  A  family  of  approximate  optimization  problems  that  are  amenable  to  computation 
is  defined  in  terms  of  least  squares  optimization  subject  to  finite  dimensional  state  space 
constraints.  We  then  give  convergence  arguments  for  approximate  optimal  parameters  to 
best  least  squares  estimates  for  the  original  infinite  dimensional  constrained  problem.  To 
demonstrate  applicability  of  our  ideas,  we  point  out  that  finite  element  approximations  in 
Galerkin  semi-discrete  formulations  based  on  piecewise  linear  spline  elements  satisfy  all  con¬ 
ditions  of  the  approximation  framework.  Included  as  a  special  case  of  our  optimal  parameter 
convergence  theory  is  the  theory  for  convergence  of  numerical  solutions  to  forward  problems 
for  nonlinear  parabolic  distributed  parameter  systems  with  delays.  We  present  sample  nu¬ 
merical  results  to  demonstrate  convergence  properties  in  both  forward  simulation  problems 
and  in  inverse  problems  with  noisy  data. 


2  Description  and  Well-Posedness  of  TCDD  Model 
2.1  The  TCDD  Model 

In  this  section,  we  present  a  mathematical  model  (1)  that  has  been  developed  [1,  2]  to  de¬ 
scribe  pharmacokinetic  and  pharmacodynamic  properties  of  TCDD.  A  convection-dispersion 
equation  (la),  based  on  the  work  of  Roberts  and  Rowland  [3],  characterizes  the  transport 
of  blood  elements  in  the  liver  sinusoidal  (blood)  region.  Throughout  this  discussion,  the  di¬ 
mensionless  spatial  variable  x  takes  values  in  the  range  [0, 1];  x  =  0  corresponds  to  the  liver 
inlet,  while  x  =  1  corresponds  to  the  outlet.  Uptake  of  dioxin  into  the  hepatic  cells,  called 
hepatocytes,  is  assumed  to  occur  by  passive  diffusion.  The  model  includes  the  dynamics 
of  TCDD-binding  with  two  intracellular  hepatic  proteins,  the  Ah  receptor  (lc)-(ld)  and  an 
inducible  microsomal  protein,  CYP1A2  (le)-(lf).  The  induction  mechanism  is  described  in 
terms  of  the  fractional  occupancy  of  the  Ah  receptor  at  a  previous  time,  t  —  Tr ,  to  account  for 
the  many  intracellular  processes  which  must  occur  before  an  increase  in  CYP1A2  concentra¬ 
tion  is  realized.  Elimination  in  the  liver  (by  metabolism  and  biliary  clearance)  is  assumed  to 
be  a  first  order  process.  A  well-mixed,  combined  venous/arterial  blood  compartment  (lg), 
which  includes  a  loss  due  to  the  uptake  and  elimination  of  TCDD  in  the  rest  of  the  body, 
completes  the  system.  A  circulatory  lag,  rc,  accounts  for  the  time  delay  in  transport  of  blood 
elements  from  the  exit  of  the  liver  to  the  venous  measurement  location. 
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The  mathematical  system  under  consideration  is  as  follows: 


(vB  +  vD^)^-  =  ~  +  p(cUH  -  fuBcB), 


luD 
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dt 
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^(t,  1)  -  QVa±dL[t.  1)  =  .4?2(f), 

C(s,x)  =  T(x),  s  G  [— r,.,0], 
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where  C  =  [Cp,  CUH,  CAh-p,  CAh,  Cpr-fy  Cpr,  Ca]T.  In  (lh),  $  and  .4g2  are  assumed  known. 

We  remark  that  gAh  and  gPr  in  equations  (lb)-(lf)  are  saturating  nonlinearities  (see 
[2])  modified  from  the  usual  product  terms, 


9Ah{y,z)  =  k+1yz,  (2) 

9pAv,z )  =  k+2yz ,  for  y,zeR,  (3) 


arising  from  the  law  of  mass  action  in  chemical  kinetics.  Specifically,  wTe  assume  that  within  a 
certain  range  of  concentrations  the  system  behaves  according  to  the  nonlinearities  prescribed 
by  (2)  and  (3)  but  eventually  saturates;  i.e. ,  due  to  the  availability  of  binding  species, 
we  assume  the  rates  of  formation  of  Ah-TCDD  complex  and  CYP1A2-TCDD  complex  are 
bounded. 


The  summary  given  above,  while  brief,  is  included  to  provide  the  reader  with  a  general 
understanding  of  the  complex  dynamics  of  the  system  under  investigation.  The  reader  is 
referred  to  the  aforementioned  works  [1,  2,  4]  for  complete  discussions. 
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2.2  Well-Posedness 


The  mathematical  model  consists  of  a  nonlinear  system  of  partial  differential  equations  with 
delays  and  questions  of  well-poseclness  are  of  interest.  A  detailed  general  theory  of  existence, 
uniqueness,  and  continuous  dependence  for  systems  which  include  the  model  above  is  given  in 
[2,  4].  The  results,  summarized  below,  were  obtained  using  a  weak  or  variational  formulation 
of  the  problem  statement  and  are  based  on  ideas  from  the  work  of  Banks  et  al.  [5]  for  nonlinear 
hyperbolic  systems. 

The  initial-boundary- value  problem  (1)  can  be  restated  in  terms  of  a  weak  or  vari¬ 
ational  formulation  (see  [2,  4]  for  details).  In  particular,  we  consider  the  state  space  V  = 
V  x  if5  x  1  and  U  =  H6  x  K,  where  V  =  H\{ 0, 1)  and  H  =  L2{ 0,1).  We  define 

Hl(0, 1)  =  1)|  m  =0}, 

with  1 -norm 

|  (j)  |  v  =  W\  h,  for  all  (j)  G  V. 

The  inner  product  in  %  is  defined  by 

6 

,  ipk)  +  fait  for  all  </>,  ril’  G  'H. 

k=  1 

where  (•,  •)  denotes  the  usual  L2  inner  product: 

{f,g)=  [  f(Z)g(Z)<%- 

J  0 

We  multiply  the  ith  equation  in  (1)  by  a  function  ^  in  a  “suitable”  class  of  test 
functions  and  integrate  in  space  in  the  first  six  equations,  followed  by  integration  by  parts 
in  the  first  equation  (la)  only.  The  weak  or  variational  formulation  of  the  problem  (1)  can 
be  written  in  operator  theoretic  form  -  see  [2,  4]  for  precise  definitions  of  the  operators  and 
nonlinear  functions  -  as  follows:  we  seek  a  solution  y(t)  G  V  satisfying 

y{t)  +  Ay{t)+AD(y{t-Tc))+g{y{t))  +  gD{y{t-Tr))  =  F{t)  in  V*,  (4a.) 

y(s)  =  y0{s),  s  £  [— r,  0].  (4b) 

Theorem  1  Under  certain  mild  assumptions  (see  [2,  4])  on  the  problem  data,  there  ex¬ 
ists  a  weak  solution  y  of  (4)  with  y  G  L2((0,T),  V)  H  C([0,T],'H)  and  y  G  L2((0,  T),  V*). 
Furthermore,  the  solution  is  unique  and  depends  continuously  on  the  data  ( yo,  F ). 

This  result  was  obtained  by  first  establishing  well-posedness  for  the  system 

y(t)  +  Ay(t)+g(y(t))  =  F(t)  in  V% 
y{s)  =  yo, 
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(5a) 

(5b) 


and  then  showing  that  the  result  could  be  extended  to  delay  systems  of  the  form  (4)  by  the 
method  of  steps  [6]  used  in  the  study  of  delay  differential  equations.  It  was  thus  noted  that 
the  system  (1),  or  more  generally,  (4),  could  be  theoretically  and  conceptually  treated  by 
investigations  of  the  abstract  system  (5). 


3  Parameter  Estimation  and  Convergence  of  Galerkin 
Approximations 


The  mathematical  model  summarized  in  Section  2.1  contains  a  number  of  physiological,  bi¬ 
ological,  and  modeling  parameters,  including  permeability  coefficients,  rate  constants,  flow 
rates,  and  the  dispersion  coefficient.  If  this  model  is  to  be  used  for  simulation  of  the  dis¬ 
tribution  and  elimination  of  TCDD  in  animals,  values  for  these  parameters  must  be  known. 
Although  some  values  can  be  measured  explicitly,  others  must  be  estimated  from  experimen¬ 
tal  data.  Thus,  a  method  to  estimate  the  so-called  ”  unknown”  parameters  is  needed. 

Based  on  the  theoretical  development  described  in  Section  2.2,  and  since  the  model 
(1)  is  a  special  case  of  (4),  which  can  be  investigated  via  (5),  we  consider  the  following  class 
of  abstract  nonlinear  parameter  dependent  parabolic  systems  evolving  in  a  real  separable 
Hilbert  space: 


y{t)  +  A{q)y{t)  +  g{q){y{t))  =  F(t;q )  (6a) 

y(0)  =  Vo-  (6b) 

In  this  case  the  unbounded  operator  A,  the  nonlinear  operator  (/,  and  the  forcing  term  F 
are  all  assumed  to  be  dependent  on  some  parameter  q.  in  contrast  to  the  ’’parameter  free” 
representation  given  in  (5). 

The  results  presented  here  are  based  on  the  general  parameter  estimation  formulation 
and  results  of  Banks  and  Kunisch  [7].  In  this  general  estimation  formulation  of  the  system  in 
(5),  we  consider  (6)  where  the  linear  operator  A.  the  nonlinear  term  g ,  and  the  input  F  have 
all  been  parameterized  by  a  vector  parameter  q  that  must  be  estimated  from  experimental 
data.  In  this  case  q  takes  on  values  from  an  admissible  parameter  set  Q.  which  may  be 
an  infinite  dimensional  set.  Given  a  data  set  of  observations  w  =  corresponding 

to  measurements  taken  at  time  tiy  we  seek  to  solve  the  general  least  squares  parameter 
estimation  problem 

I< 

min  J(q,  w)  =  V  \Cy{th  -\q)~  wd\ 2, 

<?eQ  z—' 

2=1 

where  {y{ti,  ■ ;  q)}  are  the  parameter  dependent  solutions  of  (6)  and  |  -  |  is  an  appropriately 
chosen  Euclidean  norm.  The  observation  operator  C  may  take  many  different  forms  (see  [7]) 


4 


depending  on  the  type  of  data  collected.  For  example,  if  W{  is  the  concentration  at  a  point 
z  and  time  i, .  then  C  involves  function  evaluation  in  space. 

This  minimization  problem  involves  an  infinite  dimensional  state  space  'H  and  an 
admissible  parameter  set  Q.  The  parameter  space  is  generally  infinite  dimensional  as  certain 
of  the  unknown  parameters  (for  example,  the  exit  flux  q2{t)  discussed  in  Section  2)  involve 
spatial  and/or  temporal  dependence.  We  proceed  as  in  the  works  of  Banks  et  al.  [7,  8,  9].  Let 
/Hn  be  a  finite  dimensional  subspace  of  H  and  QM  be  a  sequence  of  finite  dimensional  sets 
approximating  Q.  We  then  can  formulate  a  family  of  finite  dimensional  estimation  problems, 
with  finite  dimensional  state  spaces  and  finite  dimensional  parameter  sets,  as  follows:  find 
q  G  QM  which  minimizes 

JN(q,w)  =  | Cy*{ti,--,q)  -  v>i\ 2,  (7) 

where  yN  (t;  q)  €  'HN  is  the  solution  to  the  finite  dimensional  approximation  of  (6)  given  by: 

{yN,<f>)v*,v  +  {A(q)yN,  (p)v*,v  +  ( g{q){yN ),  4>)  =  ( F{t  q),  (p)v*,v  (8a) 

yN(0)  =  PNy0 ,  (8b) 

for  all  (j)  G  'HN ,  where  PA  is  the  orthogonal  projection  of  %  onto  %N . 

A  sequence  of  parameter  estimates  {qN,M}  results  from  the  solution  of  these  approx¬ 
imate  estimation  problems  (7)-(8).  The  question  we  address  in  this  section  is  when  one 
can  guarantee  that  this  sequence  converges  to  a  solution  of  the  original  infinite  dimensional 
parameter  estimation  problem.  A  general  unifying  framework  for  least-squares  minimization 
problems  for  first  and  second  order  systems  has  been  presented  in  [10].  Conditions  guaran¬ 
teeing  convergence  for  hyperbolic  problems  have  been  presented  in  [9]  for  linear  systems  and 
were  extended  in  [8]  to  treat  a  nonlinear  case.  We  adapt  these  results  to  a  certain  class  of 
nonlinear  parabolic  systems  which  includes  the  TCDD  model. 

In  Section  2.2  we  outlined  general  well-posedness  results  which  confirmed  conditions 
under  which  the  system  (5)  (or  (6)  without  considering  the  parameter  dependence)  has 
a  unique  weak  solution.  In  this  section  we  state  precise  conditions  under  which  (6),  as 
well  as  the  finite  dimensional  problem  (8),  have  a  unique  weak  solution  for  each  q  G  Q. 
We  also  detail  the  assumptions  on  the  general  parameter  identification  problem  which  are 
used  in  the  convergence  result,  and  state  an  abstract  sufficient  condition,  given  in  [9],  for 
the  convergence  of  the  sequence  of  solutions  {qN,M}  of  the  finite  dimensional  estimation 
problems  to  a  solution  of  the  original  infinite  dimensional  problem.  Finally,  we  show  that 
these  conditions  are  satisfied  for  a  general  class  of  nonlinear  parabolic  systems  which  includes 
the  TCDD  model. 
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3.1  Formulation  of  the  Problem 


We  assume  the  general  state  space  setting  detailed  in  [2];  that  is,  there  is  a  sequence  of  real 
separable  Hilbert  spaces  V,  "H,  V*  forming  a  Gelfand  triple  [11]  satisfying 

V  ->•  U  ~  U*  ^  V*, 

where  we  assume  that  the  embedding  V  T-L  is  dense  and  continuous  with 

\4\u  <  *#|v  for  <f)  e  V.  (9) 

We  denote  by  (•,  -)v*,v  the  usual  duality  product  [9,  11],  which  is  the  extension  by 
continuity  of  the  inner  product  in  "H,  denoted  (•,  •).  The  norm  in  %  will  be  denoted  |  •  |.  The 
operator  A(q)  is  defined  (under  the  assumptions  below)  in  terms  of  an  associated  sesquilinear 
form  a(q)  :  V  x  V  4  1;  that  is,  A{q)  E  £(V,  V*)  and  {A(q)(j),  '0)v*,v  =  c(g)(</>,  t/0- 

We  make  the  following  standing  assumptions  to  establish  our  parameter  estimation 
convergence  results.  It  should  be  noted  that  these  assumptions  are  the  same  as  those  pre¬ 
sented  in  [2]  (however  in  [2]  they  are  denoted  by  Al')-A6')),  except  that  we  now  require 
them  to  be  satisfied  uniformly  for  all  q  E  Q.  Thus,  well-posedness  of  a  weak  solution  is 
guaranteed,  under  the  conditions  of  Theorem  1  stated  in  Section  2  and  Theorem  2.2.1  in  [2], 
for  all  q  E  Q.  Our  standing  assumptions  are: 


Al)  The  form  a(q)  is  V-bounded:  for  all  <p,ip  E  V  there  exists  a  positive  constant  71 
(independent  of  q)  such  that 

M9)(<M’) I  <  7i|#]v#|v  (10) 

for  every  q  E  Q. 

A2)  The  form  a(q)  is  strictly  coercive  on  V:  for  all  (j)  E  V  there  exists  a  positive  constant 
k\  (independent  of  q)  such  that 

a(q)(o.o)  >  a  I  (11) 

for  all  q  E  Q. 

A3)  The  forcing  term  F(-\q )  satisfies 

F{-‘,  q)  E  T2((0,  T),  V*)  (12) 


for  every  q  E  Q. 
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A4)  The  function  g(q)  is  a  continuous  nonlinear  mapping  from  Li  into  Li  satisfying 

\g{q){(f>)\  <  u\d\,  (pen,  (13) 

for  some  positive  constant  7  and  for  every  q  E  Q. 

A5)  For  any  6-  A  E  V, 

{g{q){(p)  -  g(q)M,  <P  -  d)  +  hk^ld  -  A|2  >  o,  (14) 

where  k  and  k\  are  the  constants  in  (9)  and  (11). 

A6)  For  any  den  the  Frechet  derivative  of  g  exists  and  satisfies 

g’{q){(p)  e  C{Li,Li)  with  \g'{q){(p)\c(n,n)  <  C3  (15) 

for  every  q  G  Q. 

The  proofs  of  the  following  theorems  can  be  found  in  [2]  and  [4]. 

Theorem  2  Under  conditions  Al)-A6)  the  system  (6)  has  a  unique  weak  solution  y  G 
L2((0,T),Li)  for  every  yo  G  n.  The  weak  solution  depends  continuously  on  the  initial  data 
and  satisfies 

( y(t ),  <f>)vw  +  °(q){y{t),  <f>)  +  (g{q)(y(t)),  <f>)  =  ( F{t ;  g),  d)v*,v  (16) 

for  all  d  e  V  and  almost  all  t  G  (0,  T).  Moreover,  y  G  ^((O,  T),  V*). 


For  the  TCDD  model,  the  convexity  condition  A5)  011  the  nonlinearity  g(q)  is  not 
satisfied.  Weakening  this  condition  on  g,  we  established  well-poseclness  for  the  TCDD  model 
by  requiring  greater  regularity  on  the  forcing  term  F. 


Theorem  3  Suppose  V  embeds  in  'LL  compactly.  Let  F\(- ;  q)  E  L2((0,  T),Li)  and  F2(- ;  q)  G 
i/1((0,  T),  V*).  Under  assumptions  Al),  A2),  Af),  and  A6)  with  y$  E  V  and  F  =  F\  + 
F2,  there  exists  a  weak  solution  y  of  (6)  with  y  E  L2((0,T),V)  fl  C([0,T],'H)  and  y  E 
L2((0,  T),  V*).  Furthermore,  the  solution  is  unique  and  depends  continuously  on  the  data 
(yo,F). 
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3.2  The  General  Parameter  Estimation  Problem 


We  turn  to  the  least  squares  minimization  problem  described  previously:  find  q  G  Q  which 
minimizes 

K 

•/('/:  "'I  =  O/l /*•-*'/)  ~  wi\2f  (17) 

2=1 

where  are  the  parameter  dependent  solutions  of  (6)  evaluated  at  time  t,  and 

w  =  {wi}^  corresponds  to  measurements  taken  at  time  The  operator  C  is  a  continuous 
linear  operator  from  its  domain  T(J.'H)  to  the  observation  space  W,  where  J  is  a  Borel 
measurable  subset  of  the  maximal  interval  of  observations  (0,  T].  For  example,  if  we  have 
discrete-discrete  observations 

w  =  {wfaxj)}^ ti  E\V  =  Rrxl , 

then 

Cy'j  = 

where  the  natural  choice  for  the  function  space  T{J,  Fi)  is  C(J,C(  0,1)).  We  once  again 
consider  Galerkin  type  approximations  to  (6)  and,  as  discussed  in  the  introduction  to  this 
section,  define  a  family  of  approximating  parameter  estimation  problems. 

The  approximate  parameter  identification  problems  we  consider  entail  minimization 
of  (7).  As  in  [8,  9],  we  make  the  following  assumptions  for  the  infinite  dimensional  spaces 
FL,  Q  and  for  the  finite  dimensional  subspaces  FL N ,  QM: 

Bl)  The  sets  Q  and  QM  lie  in  a  metric  space  Q  with  metric  cl.  The  subspaces  Q  and 
QM  are  compact  in  this  metric  and  there  exists  a  mapping  iM  :  Q  — >■  QM  such  that 
qm  _  iM(Qy  Furthermore,  for  each  q  G  Q,  iM(q )  — >  q  in  Q  with  the  convergence 
uniform  in  q  G  Q. 

B2)  The  finite  dimensional  subspaces  l~ih  satisfy  C  V. 

B3)  For  each  i/>  G  V,  |A  —  PNip |y  — >  0  as  N  — )■  oo. 

We  further  assume  that  the  operators  A(q)  and  g(q)  as  well  as  the  forcing  term  F  depend 
continuously  on  q  G  Q;  that  is,  we  assume  that  the  following  conditions  are  satisfied: 

Cl)  |cr(g)(0t$)  -  a{q){(j),  ip)\  <  d^q,  q)  |^|v|#|v,  for  every  <f>,  F  G  V,  where  ch{q,q)  ->  0  as 

d{q,q)  ->■  0. 

C2)  | g{q){(t>)  -  g{q){F)\  <  ck{q,  q)\(/>\,  for  all  cj>  G  FI,  where  cl2{q,q)  — >  0  as  d(q,q)  — >  0. 

C3)  The  mapping  q  — )■  F(-;  q)  is  continuous  from  Q  to  L2(( 0,  T),  V*). 
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Under  the  conditions  of  Theorem  2  or  Theorem  3  and  Cl)-C3),  the  theory  of  parame¬ 
ter  dependence  of  solutions  to  ordinary  differential  equations  [12,  Theorem  6.3.1]  guarantees 
that  the  mapping  q  — >  yN(t;q )  is  continuous  for  each  t.  Furthermore,  under  condition  Bl) 
we  know  that  a  solution  {qN,M}  to  the  approximate  parameter  identification  problem  (7)-(8) 
exists.  A  general  sufficient  condition  for  the  convergence  of  {qN,M}  to  q,  a  solution  to  the 
infinite  dimensional  minimization  problem,  is  given  in  Theorem  5.1  of  [9]: 

Theorem  4  To  obtain  convergence  of  at  least  a  subsequence  of  {qN,M}  to  a  solution  q  min¬ 
imizing  (17)  subject  to  (6),  it  suffices,  under  assumption  Bl),  to  argue  that  for  arbitrary 
sequences  {qN,M}  in  QM  with  qN'M  — >■  q  g  Q,  we  have 

CyN{t]  qN,M)  — >■  Cy(t;  q).  (18) 

Note  that  condition  (18)  involves  convergence  of  the  original  sequence  of  Galerkin  approxi¬ 
mations,  which  is  a  stronger  result  than  the  subsequential  convergence  obtained  in  Chapter 
2  of  [2]  and  a  computationally  important  distinction. 


3.3  Convergence  Results 

We  show  in  this  section  that  the  convergence  criterion  of  Theorem  4  holds  under  the  general 
conditions  stated  above.  Without  loss  of  generality,  we  suppress  the  double  index  notation 
on  the  parameter  qN,M .  This  is  possible  due  to  the  assumptions  in  Bl)  -  see  [7,  9]. 


Theorem  5  Assume  that  the  conditions  of  Theorem  3  are  satisfied,  along  with  Bl)-B3)  and 
Cl)-C3).  Let  cfi  be  any  sequence  in  QN  such  that  qN  — >■  q  E  Q.  Then 


uniformly  on  [0,T].  and 


yN(t',qN )  ->■  y{t;q )  in  TL 


VN{t;qN )  ->■  y(t;q)  in  V 
for  almost  all  t.  >  0,  where  yN  satisfies 

{yN{t),  <f>)v*,v  +  o{yN){yN  (t),  <t>)  +  (. g{qN){yN  (t)),  f>) 

yN(  0)  =  PNy0 , 


for  all  <f>  G  TLN ,  and  y  satisfi.es 


{fit),  0)v*,v  +  cr(g)(y(f),  <j>)  +  ( g{q){y{t )),  4>) 

2/(0)  =  2/o, 


{F{t]  qN),  <t>)v*,v  (19a) 

(19b) 


{F(t,\q),<p)v*y  (20a) 

(20b) 


for  all  <p  G  V. 
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Proof:  We  note,  by  the  results  of  Theorem  3,  that  y(t )  G  Li  for  all  t  >  0  and  y(t)  G  V 
almost  everywhere.  By  the  triangle  inequality,  we  obtain 

I v"(P,qN)  ~  V(t-,q) |  <  | yN{t-qN)  -  PNy(t;q)\  +  \PNy(t;q)  -  y(tq)\  (21) 

for  each  t  G  [0,  T\.  Assumption  B3)  states  that  for  each  -w  G  V,  \4>  —  PA  </> |v  — >■  0  as  N  — >■  oo, 
which  of  course  implies  | A  —  PA,0|%  — >  0  for  <f>  E  Li.  Since  {y(t)  :  t  G  [0,T]}  is  compact  in 
Li.,  the  convergence 

I PNy{t-,q)  -  y{tq)\n  0 

is  uniform  on  [0,  T].  Hence,  we  need  only  show 

I yN{t]  qN )  -  -PAr2/(^;  q)\n->0  as  N  ^  oo 

uniformly  on  [0,  T]  to  obtain  the  desired  convergence  in  Li.  The  desired  convergence  in  V 
follows  in  the  same  manner  from  (21)  by  arguing 


I yN  (*,  qN )  -  pn y{t ,  9)  |v  >  0  as  n  ->■  00 

for  almost  every  t  >  0. 

Proceeding  as  in  [8],  we  define 


yN 

=  //V(/:'/Vb 

(22) 

y 

=  y{t-,q). 

(23) 

an 

=  yN(t-qN)-PNy(t-q)  =  yN  -PNy, 

(24) 

from  which  it  follows  that 

A,»  =  y?  -  lt  P»y 

(25) 

=  y?-PNyt. 

(26) 

Note  for  all  <p  G  PA , 

(Af,  0)v*,v  =  (yf  -  -^PNyA)v*y 

=  {vt  -  vt  +  yt  -  -PNy ,  4>)v*,v 

=  (yt  5  4>)v,V  —  (?/o  4>)v*,V  +  (?/*  — 

Now  by  (19a)  and  (20a)  we  have, 

(Af,0)v-,v  =  (F(tqN),(p)v,,v-a(qN)(yN^)-(g(qN)(yN),(p) 
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-  (F{t;  q),  0)y.,y  +  cr(g)(y,  <j>)  +  (g{q){y),  <j>) 

+  (?/t  — 

=  (yt  -  jtPNy,  4>')v*,v  +  (F{t;  qN )  -  F(t;  q),  0)v.,v 

-  °{qN){yN ,  0)  +  cr(g)(y,  o)  +  a{qN){y  -  PNy ,  o) 

-  ^w)(y  -  J^y, +  <y(?)(y)  -  y(?*0(yw),  0), 

=  <y*  -  ^Pwy,  </>)v*,v  +  <P(t;  ?A )  -  F(t-  q),  0)v.,v 

-  °{qN){yN  ~  PNV ,  0)  +  cr(g)(y,  0)  -  a{qN){y ,  0) 

+  <7(9")  (y  -  P"y,  0)  +  <9(9) (y)  -  yU/vS(//v).  0), 


by  adding  and  subtracting  cr(qN)(y  —  PNy,  0).  Rearranging  terms  and  using  the  definition 
of  AA  we  obtain  the  equation 


(Af,0)v%v  +  a(r¥)(AA'^) 


(yt  -  0)v*,v  +  (P(i;  )  -  F(t\  q ),  </>)v*,v 

+  cr(g)(y,  0)  -  (j{qN){y ,  (f>)  +  v{qN){y  ~  PNy ,  0) 
+  (y(y)(y)-y(yAr)(yAr),0), 


which  holds  for  all  <f>  €  P  A  .  Now,  choosing  0  =  AA  G  PA  and  recalling  that  (A^,  AA)v*,v  = 
“|AA  |2,  we  have 

~|AT  +  -(^)(A",  A")  =  <yt  -  ±P"y,  AN)V%V  (27) 

+  <F(t;  t/Ar)  -  F(t;  9),  AAV,v  +  a(9)(y,  AAf) 

—  a{qN){yi  AN)  +  cr(qN)(y  —  PNy,  AN) 

+  (g(q)(y)-g(qN)(yN),AN). 


Let  L(t)  represent  the  left  side  of  equation  (27), 

i(«)=~|AA'(()|2+ff(G)(A'Y(t),A'y)), 
and  R(t)  represent  the  right  side, 

R(t)  =  ( yt{t )  -  PNy{t ),  AN (i))v*,y  +  (P(*;  '/' )  -  P(*;  9),  AA  (t))v*,v 

+  a(q)(y(t),AN(t))  -  a(qN)(y(t),  AN(t))  +  a(qN)(y(t)  -  PNy(t),  AN(t)) 
+  (g(q)(y(t))-g(qN)(yNmAN(t)), 


which  hold  almost  everywhere.  Integrating  L  from  0  to  t  and  using  A2)  and  the  initial 
condition 


AAr(0)  =  yN( 0)  -  PNy( 0)  =  ^(0)  -  PNy0  =  0, 
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we  obtain  the  inequality 


> 


J\~\AN(s)\2  +  a(lN)(±N(s),A’'mds 


i|AK(t)|2  +  *:1  f 


AN  ( s )  |  yds 


for  alH  >  0. 


We  next  integrate  R  from  0  to  t: 


R(s)ds 

) 

=  f{(vM  -  j;pNV(s),  A s(s))v.,v  +  (F(s-,  qN)  -  F(s ;  q).  A *  (S))v.,v 

+  .(,)(■,(*),  A«W)  -  <7(,")(l/M.  A*W)  F^N)(P(»)  -  P"y(s),  AN(s)) 
+  <s(«)(i/M)  -  «(?")(!/"(»)),  Aw(*))}<fe. 


(28) 


The  first  term  under  the  integral  on  the  right  side  of  equation  (28)  is  identically  zero.  To 
see  this,  note  that  for  any  (j)  G  RN , 

(yt- ■^PNy,<P)v,v  =  (jt(y-PNy)A)v,v 

=  i<y-P«y.^, 

= 

since  y  —  PNy  G  R.  Recall  that  PN  is  the  orthogonal  projection  of  R  onto  RA ,  so  that 
(I  —  PN)y  is  orthogonal  to  elements  in  RN .  Therefore,  j-({{y  —  PN,y),(/))n  =  0  for  any 

cj)  e  RN.  "  ' 


To  treat  the  second  term  under  the  integral,  we  use  the  inequality  2 cib  <  a2  +  b 2: 

J  (F(s;q")  -  F(s\q),AN(s))v.'Vds  <  j-  J  \F(s:  qh  )  -  F(s;  q)\2v,ds 

+  e-J‘\AN(s)\ldS,  (29) 


which  holds  for  any  e  >  0. 


For  the  third  and  fourth  terms,  we  use  Assumption  Cl)  and  the  same  technique  as 
in  the  previous  argument  to  obtain 

/  {^(9)(?/(-5),  AA  (s))  -  a(qN)(y(s),  AA  (s))}ds 
Jo 
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<  [  W(q){y{s),AN{s))-a(qN){y{s)1AN{s))\ds 

Jo 

<  I  rMg,gW)|y(s)|v|Aw(s)|vds 

Jo 

<  f‘ \v(‘)\l  +  l  j‘ \&N(‘)\Us-  (30) 

Now,  using  the  assumption  that  a(q)  is  uniformly  V-bounded  for  all  q  G  Q,  we  obtain 
an  upper  bound  on  the  fifth  term  by  arguing 

[ta(qN)(y(s)-PNy(s)^AN(s))ds 

Jo 

<  71  [  \y(s)  -  P*y{s)\v\AN(s)\vds 

Jo 

<  \{]  b(s)  -  PNy{s)\lds  +  |AA  (.s)|yds.  (31) 


Finally,  we  must  obtain  a  bound  on  the  nonlinear  terms.  We  find 

[\g(q)(y(s))-g(qN)(yN(s)),AN(s))ds 

Jo 

<  r  \(g(q)(y(s))  -  g(qN)(y(s))+g(qN)(y(s))  -  g(qN)(yN(s )),  AN(s))\ds 

Jo 

<  f  \m(y(s))  -  g(qN)(y(s )),  A»)|ds  (32) 

Jo 

+  r  |<ff(?")(f/(s))  -  0(9*)(I/*(S)),  A»)|<fe. 

Jo 

To  treat  the  first  term  on  the  right  side  of  (32),  we  use  the  Cauchy-Schwarz  inequality,  (9), 
and  condition  C2),  which  states  that  g  depends  continuously  on  the  parameter  q ,  to  obtain 

r^(9)(yW)-^)(yW),A»)|da 

Jo 

<  I  \g(q)(y{s))  -  g(qN){y(s))\\AN{s)\ds 


< 

< 


<k(q,q  )|y(s)||A  (s)|ds 


k'd-Aq.  q x  ) 


2e 


|'i/(.S)|yd.S 


A:2e 

T 


AA  (s)lyd.S'. 


(33) 
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The  last  term  in  inequality  (32)  can  be  estimated  by 


l(»(9JV)(yW)-i/(9JV)(yJVW),AJV(s))|ds 


N\r„.N, 


,  N  l 


JO 


\{g{q  ){y{s))  -  g{q  ){p  y{s)) 


JO 


+g(qN)(PNy(s))-g(qN)(yN(s)),A"(s))\ds 


<  /  l(»(9JV)(y(s))-‘/(9JV)(i3iVy(a)),AJV(s))|ds 

Jo 


H 


\(g(qN)(PNy(s))  -  g(qN)(yN(s)),AN(s))\dst 


Jo 


and  is  treated  as  in  Section  2.2.5  of  [2]  when  proving  uniqueness  of  solutions.  That  is, 


l(^(^)(?/(-S'))-5(^)(JPJ"l/(.s-)),A^(S))|d.S' 


Jo 


<  f  |(  f1  g(qNY(0y(s)  +  (1  -  6)PNy(s))[y(s)  -  PNy(s)]d0,  AN(s)) \ds 
Jo  Jo 

<  C3  [t\y(s)-PNy(s)\\AN(s)\ds 

Jo 

k2C!2  C1'  k2f  P 

<  W>)-P"vM?yd>  +  ^  J  |A"(o)| Ids, 


(34) 


where  we  have  used  (9)  and  condition  A6).  Similarly, 


\(g(<f)(P"v(8))-g(<f)(vN(8)),A"(8))\d8 


Jo 


<  f  l(  /'  g(qN)'(0PNy(s)  +  (1  -  0)yN(8))[PNy(8)  -  yN(s)]dd,  AN(s))\ds 
Jo  Jo 
rt 

5iV„/  \  „,JV/„MI  aJV/ 


<  C3  \P  y(s)  —  y  (s)||A  (s)|ds 


'o 


=  C3  /  \AN(s)\2ds. 
Jo 


(35) 


Combining  (29)-(31)  and  (33)-(35),  we  obtain  the  inequality 


i|Aw(i)|2  +  [S1-(3  +  2S2)^  |A"(S)|> 

<  Y  /  lf(s;  < f)  -  F(s;  fjllv-ris  +  *  jf  |(</(s)lvfis 
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,  7i  /■*,,,  niv  ,  m2  ,  ,  k2d22{q,qN)  f\  /  u2  ? 

+  ^  y  |y(s)  - -P  y(s)ivds-f  —  y  i?/(s)iv^s 

+  — — ^-  y  |y(s)  —  PNy(s) \yds  +  C3  y  | AA  (.s)  |2c/s  (36) 

=  3 x  ( / )  ;  //  f  |Aa  (s)|2f/.s. 

Jo 

<  5n(T)  +  n  I  \AN(s)\2ds. 

Jo 

Note  that  5N(T)  — >■  0  as  A  — >■  00  by  conditions  B3)  and  C3),  the  assumption  qN  — >  q.  and 
the  result  that  y  G  L2((0,T),V).  Choosing  e  >  0  so  that  /,:  1  —  (3  +  2 A:2)|  >  0,  ignoring 
the  second  term  on  the  left  side  of  the  inequality  (36),  and  applying  Gronwall’s  lemma,  we 
obtain  the  desired  result 

A  v  ( / )  2  <  dx(T)(''r:  (37) 

that  is,  AA  — )■  0  in  C([0,T],'H)  as  N  — )■  00.  Furthermore,  ignoring  the  first  term  in  the  left 
side  of  inequality  (36)  and  using  the  result  (37),  we  find 

|AA'|V  ->•  0 

in  I/2 (0,  T)  and  hence  almost  everywhere.  Taken  with  the  remarks  at  the  beginning  of  these 
arguments,  this  establishes  the  results  of  the  theorem. 

3.4  Parameter  Estimation  for  the  TCDD  Model 

For  the  model  presented  in  Section  2.1,  each  of  the  parameters  lies  in  some  compact  subset 
of  Euclidean  space,  with  the  possible  exception  of  the  time-dependent  boundary  term,  A q-> . 
Our  initial  efforts  will  involve  a  fixed  parameterization  of  Aq->  represented  in  terms  of  a  finite 
number  of  piecewise  linear  continuous  basis  elements,  so  that  the  problem  is  conceptually 
equivalent  to  the  constant  parameter  case  in  that  the  minimization  occurs  over  a  finite 
parameter  set  [7].  In  this  case,  the  parameter  estimation  problem  will  involve  the  finite 
dimensional  approximating  sets  QM .  For  typical  examples  of  admissible  parameter  sets 
involving  both  Q  and  QM ,  the  reader  is  referred  to  [7]. 

We  state  without  proof  that  conditions  Al)-A4),  A6),  and  Cl)-C3)  hold  for  the 
TCDD  model.  This  statement  follows  from  proofs  and  arguments  given  in  [2],  For  the  state 
spaces  V  and  Ai.  conditions  B2)-B3)  are  satisfied  for  the  choice  of  piecewise  linear  continuous 
basis  elements  for  the  finite  dimensional  subspace  /HN .  Therefore,  this  example  falls  into  the 
framework  of  the  theoretical  results  presented  in  this  section. 

Finally,  we  note  that  the  TCDD  model  includes  two  time  delays,  rc  and  rrf  with 
tc  <  Tr.  A  solution  to  the  parameter  identification  problem  is  guaranteed  over  each  delay 
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interval  of  length  rc,  and  is  continued  forward  in  time  by  the  method  of  steps.  For  a  complete 
discussion,  the  reader  is  referred  to  [2], 


4  Numerical  Results 


In  this  section,  we  discuss  our  numerical  results.  First,  we  outline  the  techniques  used  in 
obtaining  a  numerical  solution  for  the  initial-boundary- value  problem  (1)  and  discuss  some 
of  the  simulated  solutions.  Then  we  show  how  we  can  use  the  simulations  for  the  forward 
problem  to  solve  the  parameter  estimation  (or  inverse)  problem  described  in  Section  3. 


4.1  Forward  Problem 

Here,  we  present  an  overview  of  the  numerical  methods  used  to  obtain  our  simulations.  All 
coefficient  matrices,  functions,  vectors,  etc.  are  as  described  in  [2],  The  values  for  the  system 
parameters  are  as  given  in  [2];  in  particular,  rc  =  1  minute  is  the  circulatory  delay,  rr  =  6 
hours  is  the  induction  delay,  and  VN  =  5.0  is  the  dispersion  number,  unless  otherwise  noted. 
Our  numerical  scheme  is  based  on  the  weak  formulation  [4]  of  problem  (1)  described  in 
Section  2. 

For  ease  of  notation,  we  define 

y  =  [C'b,  CUh  ,  C'Ah—Ti  C Ah i  Cpr-T ,  C pr ,  Ca]T  =  [yi, . . . ,  yr]T . 


Finite  Element  Formulation 


Let  0  =  x0  <  X\  <  ...  <  .z;/y  1  be  a  uniform  partition  of  the  interval  [0,1]  into  N 

subintervals  of  length  h  =  1/N.  We  take  as  basis  elements  the  standard  piecewise  linear 
continuous  functions,  4>j ,  j  =  0, . . . ,  N,  defined  by 

(  X~Xh  _1  ’  xi- 1  <:  <:  •''./• 

WO')  =  S  £2ifz£,  Xj  <  x  <  r.j .  I . 

^  0,  0  <  x  <  Xj- 1  or  xj+ 1  <  a?  <  1. 

We  define  the  Galerkin  finite  element  approximations  by 


Vi(t,x)  =  <x}(t)<j>j(x), 

y?(t.x)  =  0  (\j[l  )<> j(.r). 


(38) 


for  i  =  2, . . . ,  6,  where  the  basis  elements,  <f>j,  are  as  described  above  and  yf  is  the  approxi¬ 
mation  for  y, . 
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Next  we  substitute  the  finite  element  approximations  (38)  into  the  weak  form  of  the 
equations.  Let  yN(t )  G  ^n+5(n+i)+i  |)e  S1K.|,  that 

yN(t)  =  [a1(t),a2(t) . a«(t),yr(t)]T. 

The  finite  dimensional  system  we  obtain,  in  terms  of  the  time-dependent  coefficients  of  the 
Galerkin  approximations,  is  given  by 

MV(t)  =  ANyN{t)  +  g?(yN(t))  +  TN{t)  +  ANDyN(t  -  rc)  +  g?(yN(t  -  rr)),  (39) 

where  the  matrices  A4a  and  AN  are  elements  of  r(n+5(a+1H1)x(n+s(A'+i)+i)  ancj  the  vector¬ 
valued  functions  GNP ,  TN ,  Af}.  and  Gf  are  elements  of  MA+5(A+1)+1  (see  [2]  for  details). 

The  initial  condition,  yA ,  for  the  semi-discrete  problem  (39)  is  taken  as  the  L2  pro¬ 
jection  of  the  original  initial  condition,  <L,  onto  the  finite  element  space. 


4.1.1  General  Algorithm 

In  this  section  we  present  an  overview  of  the  general  algorithm  used  in  our  forward  problem 
simulations.  To  see  the  effects  of  the  protein  induction  delay,  r,.,  we  seek  solutions  on  the 
interval  from  zero  to  twenty-four  hours.  We  have  assumed  that  no  TCDD  is  present  in  the 
system  on  the  interval  [-6,0]  hours.  By  making  this  assumption,  we  can  ignore  the  induction 
delay  term,  C/A,  over  the  first  induction  delay  interval  [0,  r,.]. 

The  semi-discrete  problem  (39)  we  consider  is  then 

MNyN(t)  =  ANyN(t)  +  g?(yN(t))  +  TN(t)  +  ANDyN(t  -  rc)  +  gNs(yN(t  -  r,.))( 40a) 
yN{s )  =  yo(s ),  s  G  [— tc,  0].  (40b) 

The  entries  in  the  coefficient  matrices  A4A/  and  AA .  consisting  of  certain  combinations  of 
inner  products  of  basis  elements  and  their  derivatives,  were  calculated  analytically  rather 
than  using  a  numerical  integration  scheme.  The  nonlinear  vector  function  was  treated 
in  a  similar  manner.  However,  the  nonlinear  vector  function  was  calculated  using  the 
numerical  integration  scheme  quad  in  MATLAB  (The  Mat  liWorks.  Inc.,  Natick,  MA).  In 
the  following  simulations,  we  take  N  =  16;  simulations  computed  with  other  values  of  N  are 
presented  in  [2]  and  suggest  convergence  for  N  =  16. 

The  “method  of  steps”  [6]  for  ordinary  differential  equations  is  used  to  computation¬ 
ally  solve  the  problem  (40)  on  successive  circulatory  delay  intervals  of  length  rc,  and  again 
on  successive  induction  delay  intervals  of  r,..  We  first  find  the  solution  over  the  first  delay 
interval  [0.  rr\,  where  there  is  no  protein  induction  and  we  can  ignore  the  induction  delay 
term,  C?A  .  Then  we  find  the  solution  over  each  subsequent  induction  delay  interval;  on  these 
intervals,  the  protein  induction  has  begun  and  we  must  include  Gf .  In  order  to  evaluate 
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yN (t  —  tc)  and  yN (t  —  rr)  at  a  particular  time  t,  we  store  the  computed  solution  through¬ 
out  the  previous  delay  intervals  and  then  interpolate  (assuming  a  variable  step  integration 
routine  is  used)  to  find  the  value  of  the  solution  at  times  t  —  rc  and  t  —  rr.  Here,  for  ease  of 
computation,  we  assume  that  the  final  time  T  is  an  integer  multiple  of  the  induction  delay 
rr  and  the  induction  delay  is  an  integer  multiple  of  the  circulatory  delay  rc.  More  details 
regarding  the  general  algorithm  are  given  in  [2], 


4.1.2  Implementation 

The  computer  code  for  the  forward  problem  was  written  in  MATLAB  version  5.1  (The 
MathWorks,  Inc.,  Natick,  MA)  and  computations  were  carried  out  on  a  Sun  Sparc  Ultra  10 
workstation.  The  MATLAB  routine  odel5s,  which  is  a  variable  order,  variable  step  method 
based  on  numerical  differentiation  formulas,  was  used  for  time  stepping.  The  relative  and 
absolute  error  tolerances  were  set  to  1  x  10-6.  Since  this  is  a  variable  step  method,  at  time 
t  the  solution  over  the  previous  delay  interval  had  to  be  interpolated  in  order  to  determine 
the  value  of  the  solution  at  times  t  —  rc  and  t  —  rr.  In  order  to  find  the  value  of  the  solution 
at  time  t  —  rc,  MATLAB’s  interpolation  routine  interpl  was  used  with  the  spline  option. 
This  method  fits  the  data  with  cubic  splines  between  data  points,  so  that  each  segment  of 
the  curve  fit  has  at  least  the  same  first  and  second  derivatives  as  the  ones  adjoining  it.  The 
maximum  order  of  the  integrator  was  set  to  two  in  order  for  it  to  match  the  range  of  accuracy 
of  both  the  interpolation  routine  and  the  finite  element  approximations.  In  order  to  find  the 
value  of  the  solution  at  time  t  —  r,.,  linear  interpolation  was  used. 


4.1.3  Model  Simulations 


First  we  focus  on  the  solutions  on  the  time  interval  from  zero  to  twelve  hours,  in  order 
to  see  the  effect  of  the  delayed  induction  of  CYP1A2.  In  these  simulations  we  consider  a 
special  case  of  the  original  model  to  facilitate  the  computation  of  simulated  solutions  and  to 
make  the  effects  of  the  induction  delay  term  more  pronounced.  We  use  a  special  case  of  the 
boundary  condition  at  x  =  1;  we  set  v  =  0  and  Q2{t)  =  0  in  our  original  boundary  condition 

BC 

vCB(t,  1)  +  Ca(t)  -  1))  =  q2(t). 


to  obtain 


BCp 

dx 


(M)  =  0. 


This  condition  permits  computation  of  simulated  solutions  over  longer  time  intervals  without 
error  accumulation  so  that  qualitative  properties  can  be  demonstrated.  Since  it  is  a  special 
case  of  our  original  boundary  condition,  the  theory  discussed  previously  still  holds.  In 
addition,  we  set  IPr,  the  maximum  induction  rate  of  CYP1A2,  equal  to  434.6  nmol/L/hr 
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which  is  ten  times  the  rate  suggested  by  the  literature  [13].  This  is  done  to  magnify  the  effects 
of  the  induction  delay  in  our  computed  solutions.  Preliminary  computational  investigations 
on  the  qualitative  behavior  of  the  system  without  these  changes  are  discussed  in  [14]  and  [2], 
In  particular,  the  behavior  of  the  system  on  the  time  interval  from  zero  to  six  hours  and  the 
dependence  of  the  solution  on  the  dispersion  number,  VN ,  and  the  circulatory  delay,  rc,  are 
presented. 

During  the  time  interval  from  zero  to  six  hours  (when  t  <  Tr).  the  protein  CYP1A2 
is  present  only  at  a  basal  level.  After  six  hours  (t.  =  rr)  however,  there  is  a  dose-dependent 
induction  of  CYP1A2.  This  induction  is  represented  in  the  model  by  the  induction  term 
C/f  {y{t  —  r,,)).  Since  the  model  incorporates  the  binding  of  CYP1A2  with  TCDD,  as  well 
as  the  other  dynamics  of  TCDD  interaction  in  the  liver,  we  expect  that  the  inclusion  of  the 
induction  term  in  the  model  affects  in  a  nonnegligible  way  the  solutions  of  the  system.  We 
compare  the  simulated  solutions  with  and  without  the  induction  term  to  study  the  nature 
of  these  effects. 

First  we  consider  the  concentration  of  available  CYP1A2  in  the  hepatocytes,  denoted 
Cpr  (Figure  1).  We  note  a  sharp  increase,  followed  by  a  gradual  decline,  in  C[p,.due  to  the 
beginning  of  the  induction  of  CYP1A2  and  its  subsequent  binding  with  TCDD.  The  effect  of 
this  binding  can  also  be  seen  in  Cpr-p ,  the  concentration  of  CYP1A2-TCDD  complex  in  the 
hepatocytes  (Figure  2).  Here  we  note  an  initial  period  of  adjustment  to  the  elevated  levels 
of  CYP1A2,  and  then  an  increase  in  Cpr-p  due  to  the  binding  of  CYP1A2  with  TCDD.  In 
addition,  the  binding  of  TCDD  results  in  a  relative  decrease  in  the  concentration  of  unbound 
TCDD  in  the  hepatocytes,  CUH  (Figure  3).  Since  there  is  a  decrease  in  the  concentration  of 
unbound  TCDD  in  the  hepatocytes,  unbound  TCDD  in  the  liver  blood  diffuses  into  the 
hepatocytes,  where  the  CYP1A2  is  produced;  this  results  in  a  relative  decrease  of  both  CUB , 
the  concentration  of  unbound  TCDD  in  the  liver  blood,  and  Cp,  the  total  concentration  of 
TCDD  in  the  liver  blood  (Figures  4  and  5).  Furthermore,  the  induction  of  CYP1A2  indirectly 
affects  the  concentration  of  the  Ah  receptor,  the  other  protein  with  which  dioxin  binds.  We 
note  that  the  relative  decrease  in  CUH  causes  a  relative  increase  in  C^,  the  concentration 
of  available  Ah  receptor  protein  in  the  hepatocytes,  and  a  relative  decrease  in  C,\h-T ,  the 
concentration  of  Ah  receptor-TCDD  complex  in  the  hepatocytes  (Figures  6  and  7). 


4.2  Inverse  Problem 

As  we  mentioned  at  the  beginning  of  Section  3,  our  model  contains  many  physiological, 
biological,  and  modelling  parameters.  When  we  computed  the  forward  problem  simulations, 
we  set  these  parameters  using  values  given  in  the  literature.  In  order  to  use  this  model 
to  simulate  the  disposition  and  elimination  of  TCDD  in  a  specific  animal,  we  must  know 
the  values  for  the  parameters  for  the  animal  in  question.  Many  of  these  values  must  be 
estimated  from  experimental  data.  In  Section  3,  we  discussed  the  theory  behind  such  a 
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parameter  estimation  problem;  here  we  provide  numerical  results  for  a  problem  of  this  type. 
For  a  demonstration  example,  we  estimate  the  value  of  the  axial  dispersion  number  VN.  a 
parameter  that  is  not  directly  measurable  and  hence  one  that  will  not  generally  be  available 
from  experimental  literature.  We  do  this  holding  all  other  parameter  values  fixed.  Since 
we  do  not  have  access  to  experimental  data,  we  use  the  solution  for  the  arterial  blood 
concentration  Ca  from  the  forward  problem  simulations  as  our  data.  This  means  that,  in  the 
context  of  Section  3,  the  observation  operator  C  is  the  dot  product  of  the  solution  at  time  t 
with  the  unit  row  vector  in  Rlx(A+5(A'+i)+i)  with  a  one  in  the  last  component. 


4.2.1  General  Algorithm 

In  this  section,  we  define  our  notation  and  outline  the  general  algorithm  used  in  our  simula¬ 
tions.  Our  goal  is  to  find  the  value  for  VN  that  minimizes  the  difference  (in  the  least  squares 
sense)  between  Ca(t,VN),  the  ©^-dependent  solution  from  the  forward  problem,  and  our 
data.  We  let  Ca{ti ,  VN)  be  the  simulated  solution  of  Ca  at  time  £*,  with  t,  in  the  time  interval 
from  zero  to  three  hours.  We  denote  the  vector  of  data  points  or  observations  by  Cn .  and 
we  let  ( Ca)i ,  the  value  of  Ca  observed  at  time  U,  be  the  zth  component  of  Ca.  Finally,  C*  is 
the  simulated  solution  of  Ca with  VN  =  V*N  =  5.0;  i.e.,  we  take  V*N  =  5.0  as  the  true  value 
in  our  simulated  data. 

The  general  algorithm  consists  of  defining  a  least  squares  cost  function  and  choosing 
an  initial  value  (or  values)  for  VN ,  the  parameter  to  be  estimated,  which  are  then  used  in 
an  optimization  routine.  Here  we  use  a  Nelcler  Mead  optimization  routine  [15].  Since  Nelder 
Mead  is  a  simplex  search  algorithm,  we  are  required  to  choose  an  initial  simplex  for  DiV, 
which  in  our  case  can  be  thought  of  as  an  ordered  pair  of  initial  values. 


4.2.2  Implementation  and  Data  Creation 

The  algorithm  and  inverse  problem  feasibility  was  tested  with  both  noise  free  data  and 
data  containing  various  levels  of  random  noise.  The  necessary  computer  code  was  written 
in  MATLAB  version  5.1,  and  the  computations  were  carried  out  on  a  Sun  Sparc  Ultra  10 
workstation.  In  the  Nelcler-Mead  simplex  search  algorithm  for  the  optimization  [15],  we  set 
the  termination  tolerance  at  10-4  when  no  noise  was  present  in  the  data  and  at  10-3  in  the 
presence  of  noise. 

In  a  lab  setting,  we  would  expect  to  see  noisy  data  due  to  such  factors  as  measurement 
error,  human  error,  and  environmental  effects.  Since  we  did  not  have  experimental  data  for 
Ca,  we  use  three  types  of  simulated  data  in  our  inverse  problem:  data  without  noise,  data 
with  uniformly  distributed  random  noise,  and  data  with  normally  distributed  random  noise. 
We  compared  our  ability  to  estimate  the  value  of  the  parameter  when  we  added  1%  relative 
error  and  5%  relative  error  so  we  could  demonstrate  the  effect  of  the  magnitude  of  the  noise 
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on  our  ability  to  estimate  the  parameter  correctly.  All  three  data  sets  were  derived  from 
C*(t ),  the  simulated  solution  with  the  “correct”  value  VN  =  V*N  =  5.0.  Below  we  describe 
how  the  data  sets  were  created. 

The  data  without  noise  is  simply  the  original  time  domain  solution  C*(t)  (Figure  8). 
We  added  uniformly  distributed  random  noise  to  the  original  time  domain  solution  so  that 
the  “data”  Cau  (Figures  9  and  10)  was  given  by  ( ( „ ) ,•  =  (1  +  ( -,)( ',](/,)•  In  this  case,  the 
e,;  were  uniformly  distributed  random  numbers  that  were  appropriately  scaled  and  shifted, 
depending  on  the  magnitude  of  the  noise,  so  that  c,  G  [—.01,. 01]  for  1%  relative  error  and 
c,  G  [—.05,  .05]  for  5%  relative  error.  The  random  numbers  were  generated  by  the  MATLAB 
function  rand ,  which  outputs  uniformly  distributed  random  numbers  from  0  to  1.  We  added 
normally  distributed  random  noise  to  the  original  solution  to  create  Can  (Figures  11  and  12) 
in  a  similiar  fashion,  so  that  ( Can)i  =  (1  +  ej)C*(£j).  Here,  the  e,;  were  normally  distributed 
random  numbers  with  mean  0  and  variance  1  that  were  scaled  accordingly,  so  that  with  99.7% 
certainty  e*  G  [—.01,  .01]  for  1%  relative  error  and  e*  G  [—.05,  .05]  for  5%  relative  error.  The 
MATLAB  function  randn  generated  these  normally  distributed  random  numbers. 


4.2.3  Model  Simulations 

Using  the  data  and  methods  described  previously,  we  performed  numerous  inverse  problem 
calculations  to  investigate  our  ability  to  estimate  the  value  of  the  parameter  VN.  We  con¬ 
sidered  several  different  situations:  using  “correct”  data  and  a  range  of  initial  conditions, 
using  the  data  with  uniformly  distributed  random  noise  adding  1%  and  5%  relative  error, 
and  using  the  data  with  normally  distributed  random  noise  adding  1%  and  5%  relative  error. 

We  find  that  when  we  use  the  data  set  without  noise,  we  are  able  to  recover  the  true 
value  of  VN  ( V*N  =  5.0)  with  at  least  four  decimal  places  of  accuracy  and  obtain  cost  function 
values  on  the  order  of  10-5.  This  is  true  for  initial  estimates  for  VN  with  components  in  the 
range  [4.5,  5.5].  Since  there  is  no  noise  in  the  data,  these  results  are  what  we  would  expect 
and  completely  typical  of  many  inverse  problem  calculations  using  data  without  noise.  In 
Table  1  we  present  typical  findings. 


Initial  Simplex 

x0 

Converged  Value 
VN 

True  Value 
V* 

UN 

Cost  Function 

J 

(4.5,4.51) 

5.0000 

5.0 

4.7301  x  10“5 

(5.5,5.51) 

5.0000 

5.0 

4.7301  x  10“5 

Table  1:  Results  for  data  without  noise. 


When  we  use  the  data  with  noise,  we  find  that  we  are  still  able  to  solve  the  inverse 
problem  successfully.  Even  with  up  to  fifty  percent  error  in  our  initial  estimate,  the  converged 
values  VN  obtained  from  the  optimization  routine  are,  at  worst,  equal  to  V*N  with  error  on 


21 


the  order  of  10-3  and  at  best,  are  equal  V*N  with  at  least  four  decimal  places  accuracy.  The 
values  of  the  cost  function  range  from  &  153  to  ~  11,609.  Tables  2  and  3  illustrate  some  of 
these  results. 


%!  Noise 

Initial  Simplex 

,T0 

Converged  Value 
Vn 

True  Value 
T>* 

UN 

Cost  Function 

J 

1% 

(5.0,5.01) 

4.9990 

5.0 

4.5320  x  102 

1% 

(4.5,4.51) 

4.9990 

5.0 

4.5320  x  102 

1% 

(5.5,5.51) 

4.9990 

5.0 

4.7249  x  102 

5% 

(5.0,5.01) 

5.0034 

5.0 

1.1609  x  104 

5% 

(4.5,4.51) 

5.0034 

5.0 

1.1609  x  104 

5%, 

(5.5,5.51) 

5.0034 

5.0 

1.1609  x  104 

Table  2:  Results  for  data  with  uniformly  distributed  noise. 


%  Noise 

Initial  Simplex 

,T0 

Converged  Value 
Vn 

True  Value 
T)* 

UN 

Cost  Function 

J 

1%. 

(5.0,5.01) 

5.0000 

5.0 

1.5342  x  102 

1% 

(4.5,4.51) 

5.0000 

5.0 

1.5342  x  102 

1% 

(5.5,5.51) 

5.0000 

5.0 

1.5342  x  102 

1% 

(2.5,2.51) 

5.0000 

5.0 

1.5342  x  102 

1% 

(7.5,7.51) 

5.0000 

5.0 

1.5342  x  102 

5% 

(5.0,5.01) 

4.9998 

5.0 

3.8822  x  103 

5% 

(4.5,4.51) 

4.9998 

5.0 

3.8822  x  103 

5% 

(5.5,5.51) 

4.9998 

5.0 

3.8822  x  103 

Table  3:  Results  for  data  with  normally  distributed  noise. 


We  find  that  we  do  a  better  job  of  recovering  the  true  value  V*N  when  there  is  normally 
distributed  noise  than  when  there  is  uniformly  distributed  noise,  no  matter  what  the  percent 
error.  This  is  to  be  expected  since  the  noise  in  the  data  is  more  likely  to  be  “larger”  with  the 
uniformly  distributed  noise.  To  see  this  more  clearly,  we  can  look  at  e,  for  any  i,  which  we 
described  in  the  previous  section  when  we  discussed  the  creation  of  the  data.  We  consider 
the  case  with  1%  relative  error.  For  the  uniformly  distributed  noise,  it  is  equally  likely  that 
€j  is  any  number  in  [-.01,  .01].  However,  for  the  normally  distributed  noise,  there  is  a  68% 
chance  that  the  random  number  generated  by  the  routine  is  within  one  standard  deviation 
(in  our  case  a  =  1)  from  zero  and  a  99.7%!  chance  that  it  is  within  three  standard  deviations 
from  zero.  Here,  that  means  that  there  is  a  99.7%  chance  that  e,  £  [—.01,  .01]  and  a  68%) 
chance  that  e,  £  [—.0033,  .0033].  So,  it  is  likely  that  the  normally  distributed  noise  has  a 
smaller  magnitude  than  the  uniformly  distributed  noise. 
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5  Concluding  Remarks 


The  previous  discussions  provide  a  rigorous  foundation  for  computational  methods  to  in¬ 
vestigate  both  forward  simulations  and  inverse  problems  for  nonlinear  parabolic  partial  dif¬ 
ferential  systems  with  delays.  In  particular,  the  motivating  TCDD  spatially  distributed 
liver  model  can  be  treated  using  the  computational  framework  presented  here.  The  sample 
computational  results  for  this  model  demonstrate  the  applicability  of  our  approach  in  the 
investigation  of  both  qualitative  and  quantitative  properties  of  such  models. 
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Figure  1:  The  effect  of  the  induction  delay  on  the  concentration  of  available  CYP1A2  in  the 
hepatocytes. 


Figure  2:  The  effect  of  the  induction  delay  on  the  concentration  of  CYP1A2-TCDD  complex 
in  the  hepatocytes. 
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Figure  3:  The  effect  of  the  induction  delay  on  the  concentration  of  unbound  TCDD  in  the 
hepatocytes. 


Figure  4:  The  effect  of  the  induction  delay  on  the  concentration  of  unbound  TCDD  in  the 
liver  blood. 
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Figure  5:  The  effect  of  the  induction  delay  on  the  total  concentration  of  TCDD  in  the  liver 
blood. 
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Figure  6:  The  effect  of  the  induction  delay  on  the  concentration  of  available  Ah  receptor  in 
the  hepatocytes. 


Figure  7:  The  effect  of  the  induction  delay  on  the  concentration  of  Ah  receptor-TCDD 
complex  in  the  hepatocytes. 


Figure  8:  The  computed  solution  for  the  concentration  of  TCDD  in  the  arterial/venous  blood 
compartment  ( Ca )  without  noise. 
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Figure  9:  The  computed  solution  for  the  concentration  of  TCDD  in  the  arterial/venous  blood 
compartment  (Ca)  with  1%  uniformly  distributed  random  noise. 


Figure  10:  The  computed  solution  for  the  concentration  of  TCDD  in  the  arterial/venous 
blood  compartment  (Ca)  with  5%  uniformly  distributed  random  noise. 
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Figure  11:  The  computed  solution  for  the  concentration  of  TCDD  in  the  arterial/venous 
blood  compartment  (Ca)  with  1%  normally  distributed  random  noise  (99.7%  certainty). 


Figure  12:  The  computed  solution  for  the  concentration  of  TCDD  in  the  arterial/venous 
blood  compartment  (Ca)  with  5%  normally  distributed  random  noise  (99.7%  certainty). 
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